#include "MPC_C_T.h"

MPC_C_T::MPC_C_T(	double	iRough2_set,
					double	iFact_R2, 
					double	iM2_set,
					double	iM2_fact,
					double  iT_Low,
					double  iT_High,
					double	iRT_T,
					double	iT_ini,
					double  iW_default
					double	it_end,
					double	it_int_ctrl, 
					int		iStep_num,
					int		iRel_err,
					int		iDetail):
MPC_Base( 0.98,
		  0.0,
		  iRough2_set,
		  iFact_R2,
		  300,
		  0.0,
		  it_end,
		  iT_ini,
		  iRT_T,
		  it_int_ctrl,
		  iStep_num,
		  iRel_err,
		  "MPC_C_T.dat",
		  iDetail,
		  1)
{
	;
}
// ---------------------------------------------------------------------------
double MPC_C_T::update(double time,lattice_state iLState,const double P, const double T)
{
	if(time>=t_arr[t_idx-1])// Don't need an update
		return T;
	else // Solve optimization problem
	{
		
		int i;
		int ndim,nfunk,paralength = 1;
		double *temp,*para,*uplim,*lolim;
		double ftol = 1e-8;
		double temp_upper,temp_lower;
		double temp_ini;
		double cost_min,cost_temp,err;

		if(iLState.SOR>1.0 || iLState.SOR < 0.0)	{cerr << "\nError@MPC_Base::update(): illegal value of density. iLState.SOR = " << iLState.SOR << "\n";system("pause");exit(0);}
		if(iLState.R2 < 0.0)						{cerr << "\nError@MPC_Base::update(): illegal value of roughness. iLState.R2 = " << iLState.R2 << "\n";system("pause");exit(0);}
		if(iLState.H < 0.0)							{cerr << "\nError@MPC_Base::update(): illegal value of height. iLState.H = " << iLState.H << "\n";system("pause");exit(0);}
		if(P <= 0.0)								{cerr << "\nError@MPC_Base::update(): illegal value of pressure. P = " << P << "\n";system("pause");exit(0);}
		if(T <= 273)								{cerr << "\nError@MPC_Base::update(): illegal value of temperature. T = " << T << "\n";system("pause");exit(0);}		

		alphas	= iLState.alphas;
		betas	= iLState.betas;
		ht		= iLState.H;
		R2		= iLState.R2;
		g_dens	= iLState.SOR;

		para=vector(1,1);
		para[1] = 1;
		temp  = vector(1,1);
		uplim = vector(1,1);
		lolim = vector(1,1);
		ndim  = 1;

		cost_min = 10e20.0;
		for(temp[1] = T_low;temp[1] <T_high;temp[1]++)
		{
			cost_temp = cost_funct(temp,para);
			if(cost_temp<cost_min)
			{
				cost_min = cost_temp;
				temp_ini = temp[1];
			}
		}
		free_vector(temp,1,1);

		para[1] = Step_num;
		temp = vector(1,Step_num);
		uplim = vector(1,Step_num);
		lolim = vector(1,Step_num);
		ndim = Step_num;



		for(i=1;i<Step_num;i++)
		{
			temp[i] = temp_ini;
			temp_upper = T+i*t_int_ctrl*RT_T;
			temp_lower = T-i*t_int_ctrl*RT_T;
			uplim[i] = T_high>temp_upper?temp_upper:T_high;
			lolim[i] = T_low>temp_lower?T_low:temp_lower;
		}

		err = fminsearch(temp,ndim,ftol,cost_funct,para,&nfunk,uplim,lolim);
		for(i=1;i<Step_num;i++)
		{
			temp[i] = uplim[i]>temp[i]?temp[i]:uplim[i];
			temp[i] = lolim[i]>temp[i]?lolim[i]:temp[i];
		}

		T = temp[1];]
		free_vector(lolim,1,Step_num);
		free_vector(uplim,1,Step_num);
		free_vector(temp,1,Step_num);
		free_vector(para,1,1);

		t_idx ++;

		return T;
	}	
}